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SUI-31ARY 


Several tasks pertinent to flight control in parameter uncertainty 
and wind-gust loading have been successfully completed. 

Identification Algoritlms fcr extracting stability and control 
derivatives from flight data taking gust loading into account have been 
developed. They have been verified by simulation and evaluated throughly 
on actual flight data taken on a Lockheed Jet Star flying in turbulence. 

In particular the need for automatically generated dither-like inputs has 
been studied. 

Criteria for performance evaluation using stochastic models have been 
developed for gust alleviation as well as handling qualities. Algorithms 
for assessing degradation in performance due to parameter uncertainty 
have been developed and evaluated using flight test data (Lockheed Jet 
Star). 



INTRODUCTION 


This report deals with the general problem of flight control taking 
into account parameter uncertainty and wind gust loading. The following 
specific tasks were undertaken: 


Task A. Develop criteria for assessment of control 
performance for uncertain systems with par- 
ticular reference to flight control under 
parameter uncertainty and gust loads. 


Task B. Develop identification algorithms for 

stability and control derivatives from flight 
data taking gust loading into account , and 
assessing the need for automatically generated 
dither-like inputs. 

Task C. Evaluate the performance of identification 
algorithms for typical flight conditions. 


The report is organized generally in order of these tasks. 

V. T e begin in section 1 with the- problem, of control performance 
assessment. V.'e consider two kinds of performance evaluation: find- 

Alleviation and Handling Quality Improvement. Vie develop a general 
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in section 1.1 which embraces indices of performance for both categories, 
and present formulas for assessing degradation in performance due to 
parameter uncertainty. In section 1.2 we illustrate the theory by 
evaluating the degradation in gust alleviation for typical flight condi- 
tions using Lockheed Jet Star data. In particualr we use stability and 
control parameters obtained with dither-like inputs as reflecting uncer- 
tainty in parameters , and in this way assess the efficacy of dither-like 
inputs for parameter identification in flight test data containing gust 
loading. A new handling quality criterion is developed in section 1.3 
vising a stochastic pilot input model (see L. W. Taylor [ll ) as well as 
the corresponding optimal feedback control. 

Section 2 is devoted to Task B. Using a continous-time white noise 
theory developed by the Principal Investigator, a maximum likelihood 
identification algorithm for flight test data containing wind gust response 
is developed in sections 2.1 and 2.2. The first results using this theory' 
were reported by K. W. Iliff in his UCLA thesis [2]. The specifics of the 
algorithms employed in the present work are quite different however. The 
algorithms were verified by simulating longitudinal short, period motion 
accounting for vertical wind gust modelled by the Dryden Spectrum. The 
simulation results are presented in section 2.3. 

The bulk of the results presented in this report is contained in 
section 3 where we evaluate the performance of identification algorithms 
(developed in section 2.2) for typical flight conditions, as required 
under task C. Flight test data taken on a Lockheed Jet Star flying in 
heavy turbulence (obtained courtesy of DFRC) was used. The tine-histcrv 
of about 60 seconds was divided into several sections and the data in each 
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section analysed separately. The stability and control derivatives were 
extracted as well as the gust intensity. Moreover the vertical gust 
velocity wave-form has been reconstructed and the corresponding power 
density spectrum plotted. The results generally agree with those of [2], 
and provide a thorough and extensive evaluation of the identification 
algorithms for typical flight conditions. 

Section 4 utilized the Jet Star flight data to study the efficacy of 
using dither-like inputs in place of conscious pilot inputs. Two meneuvers 
not containing pilot input were used for this purpose and the algorithms 
used to determine the derivatives. The extracted values are compared with 
the derivatives obtained during an adjacent maneuver with a significant 
pilot input. 

All the tasks. A, B, and C have been successfully completed, as the 
report will indicate. 

Fi nall y, and Appendix (Appendix I) examines a basic question in 
flight -test data processing: whether system models should be time- 

discrete or time-continuous. It is shown that for type of flight-test 
data considered the continuous-time model is to be preferred. 
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1. ASSESSMENT OF CONTROL PERFORMANCE UNDER RAPAITTEP. UNCERTAHT/ AND GUST 
LOADING. 

1.1 The Basic Formulas 

In this section we develop some general formulas for the assessment 
of the degradation in control performance due to parameter uncertainty. 

As we shall show in Section 3, it is general enough in formulation to 
i ncl ude both gust alleviation and handling quality criteria. 

We begin with th- linearized state equations of motion: 


x(t) = Ax(t) + Bu(t) + F n(t) 
y(t) = Cx(t) + G n(t) 



( 1 . 1 ) 


where x(t) is the state (enhanced as necessary to take care of the wind 
gust model with the Dryden or similar rational spectrum) , and y(t) re- 
presents the sensor data available. The matrix A represents the stability 
derivatives, and the matrix B contains the control derivative;, 

G G* = I 
F G* = 0 , 

and n(t) is white Gaussian with unit spectral dans: ty . The performance- 
criterion is: 

n 

(l/T ) O 7 

Jo (i:(! |i.'-(t) 1 1" + ? Iiu(t)| | ) dt (i.:-) 

where L is a matrix that determines the variable to be minimised (e.r. t: ,? 


normal acceleration), and A is a suitably chosen positive number to yield 
the appropriate level of control constraint. E denotes expected value. 

We are interested primarily in the case where T is large (or, T , in 
theory) . 

We know that the optimal control that minimises the porfcrmar.ce 
criterion chosen is given by: 

u (t) = -B* P x (t) /A (1.3) 

O CO 

where P c satisfies the Control Riccatrti Equation: 

P (t) + A* P (t) + P (t)A + L*L - P (t)BS*P (t) /A =0 
c - c c c C 

r / ' 

P (T) = 0 
c 

As we have indicated, we are primarily interested in the case of 'large' 
T, so that we may use the ’steady state' version, setting the time 
derivative to be zero: 

A : ' : ? + P A + L*L -PE E ; ' : ? /A = 0 (l.u) 

c c c c 

The Kalman estimate x(t) is defined by: 

x(t) = Ax(t) + B u (t) +F C C*(y(T)-Cx(t) ) ( 1 . 5 ) 

C x 

where, consistent with out steady state view-pcir.t, we shall take the 
initial estimate x(0) to be zero, arid F r the filter matrix satisfies: 



A p f + p f A* +FF* - P f C* C P f = 0 


( 1 . 6 ) 


For the optimal system, the performance functional 


where 


J = Lim 
T 


ft) i: 


|L x(t) | | dt 


= Tr (L*LP- + L* L J ) 
f a 


(1.7) 


(A-BB'-'-'Pc )J + J (A-BB*P„ ) + P f C : ' : CP = 0 

1 — - A Pi ; C I I 


A a a A 


Of course, to achieve the value of the performance criterion given 
by (1.7) we will need to know the systems parameter values exactly. Our 
task now is to determine what happens when the parameters are unknown, and 
have to be estimated by an identification algorithm, or other means. Let 
us denote the estimated parameter matrices by the subscript e. Then the 
feedback control will be determined as: 


u (t) = -B *(P J A ) x(t) 
o e c,e 


0.3) 


where 


A* P + P A + L * L - (1/A) F B F * T = 0 0.7) 

e c,e c,e e e e c,e e o c,e 


x(t) = A x(t) + I? u (t) + F C * ( v ( 1 ) - C x(t) 
e e o f ,e e • c 


} ( 1 . 1=0 


X(0) 


0 . 



r 

t 


A P* + P. A * + r F * - P r C C * !’_ =0 0 

e f,e f,e e e e f,e e e f,e (1.11) 


The problem now is to calculate the functional J for thin choice of 
optimal control. We know of course flint it will be lorpcr than that p.ivcn 
by (1.7). We can calculate it as follows, let 


z(t) = 


x(t) 

A 

x(t) 


( 1 . 12 ) 


Under the choice (1.8) thru (1.11), we have: 


x(t) 

x(t) 


= A x(t) - (1/A) B B * T x(t) + F n(t) 

e c ,e 


(1.13) 


(A - P. C* C ) X(t) + P, C* ( C x(t) + G n(t) ) 
e f,e e e t,e e ■ 

- B B * P x(t)/A (1.14) 
e e c,e 


where we substituted for y(t) fixam (1.1). We can rein' it e these equations 
in the form: 


z(t) = A "(t) + H 

1 *e e 


n(t) 


where A e is the compound matrix: 


\ - 


A 


r> c*c 

f ,e e 


-BP- P A 
e c.e 


A 


F, C" C -B B"P />. 
1 ,e e e o e c,c 


0 . 1 : > 


P. C* G 
1 ,e e 


-'i- 



/ 


Then letting K denote the compound matrix: 
K = L 0 

and 

M = 0 I 

so that: 

L x(t) = K z(t) 

M z(t) = x(t) 

and defining 

R - Lim E( z(t) z(t) : ' : ) 
e 


we have 


J = lim (1/T) I | |Lx(t)| | dt 
0 


Ir. K*K R (1.17) 
e 


Of course this simplifies to (1.7) when the parameters are correct, he 
note that R g may be computed by solving the linear equation: 


A R + 
n e e 


r A“ + H H 
e 'e e e 


0 


The actual degradation in performance is of course the diflerer.ee 
between (1.17) and (1.7). The calculation however can be made more 
efficient by noting that the difference will be small and is well ap- 
proximated by the linear terms only. In other words, denoting the 
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difference by AN, we have: 

AN = Tr. (K*K AR g ) 


inhere AR can be calculated by the linear approximation equation: 
e 


A aR + (AR )A" + <aA)R + R(aA/ 


H(AH)* 


(AH)H* 


And AP is calculated fran: 
c 

(A* - (1/A)P C B B*)(AP c > + Adjoint . + L* AL + (AL) :: 

+(AA)*P + P c AA - P C (AB)B* P c /X ~ Adjoint = 0 

Similarly, AP f is obtained from the liner equation: 


(A - P f C C*)AP f + Adjoint + ( (AF)* + F (AF)- 

+(AA)P f + P f (AA)* = 0 


- P (AC) C* P f 


= 0 ( 1 . 20 ) 


L 

( 1 . 21 ) 


Adjoint 

( 1 . 22 ) 
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1.2 Application to Jet Star 

We now specialise the formulas of the previous section to calculate 
the degradation in control performance due to uncertainty in the parameters 
to the case of the Lockheed Jet Star. We consider the linearised longitu- 
dinal short period motion equations: 


Z a + 0 + z 6 + Z w 

a oe e a g 


M a + MtO + M. 6 + M w 

a 0 6e e a g 


-20 6 + 20 <5 

e p 


where 6 ^ is the elevetor deflection due to the pilot input 6 p (we are 
neglecting the servo-actuator non-linearity), and the other quantities 
are as before. We assume a Dryden spectrum for the vertical component 
of the gust velocity , and hence we can write : 


W 5 

_£ 


-v a + Oj ( /2|vL > N(t) 

h 


where N(* ) is white Gaussian with unit spectral density, is the scale 
lenght of the turbulence. The break frequency (radians |scc) is: 


w = v 


- 1 ?- 



and the intensity is: 


We take 


4 = 


= 1000 


in this section. 

The sensor measurements are: 


q a Pitch Pate 
0 a Pitch angle 


n 'u Normal acceleration 
z 

'v. a-Vane deflection angle 


we note that 


(0- o] 


K (a + w ) 
“ 

v 


K a £ a G 


VJe can put these relations in State Space form, allowing for sensor measu 
ment errors at the same time. Vfe take the errors to be white Gaussian. 


- 1 ?- 


Let 



where the subscrip m denotes the senscr measurement corrupted by noise. 
Then we can write 

v = C x + G N 
x = Ax+Bu+FN 


-1‘j- 






2 2 2 2 

where g-j_, g^, g^ are white noise spectral densities corresponding to 
the corresponding sensors. 

We take our performance criterion to be that of minimising the normal 


acceleration due to gust. Hence 



and is lx 5. 

The nominal parameter values wall be taken to be those listed below, 
corresponding to those extracted from maneuver A. 


- if. - 



[ 


Z a 

zz 

-1.45 

H 

a 

- 

-9.79 

M 

q 

- 

-1.43 

Z « 

e 

= 

-0.098 

M 6 

e 

— 

-8.15 

2 

°1 

= 

25.2 



The corresponding value of m.s . normal acceleration J was calculated fcr 
two values of A : A = 1 (large control), A = 10 (small control). Fcr 
A = 1, J = .007985; for A = 10, J = .008264. It nay be noted that in the 
absence of any feedback control, J = .01185. The reduction is not spectacu- 
lar for the particular control configuration, but our purpose here is to 
illustrate the techniques. To calculate the performance degradation (due 
to wrong parameter choice) we shall take the parameter' values extracted 
from maneuver X corresponding to the 'dither' input. In other words, wo 
are determining the adequacy of parasieters extracted using only a dither 
input (in gust, of course). Hie parameter values corresponding to maneuver 
X are. 


Z = -1.41 
a 

M = -9.42 

a 

M = -0.53 

q 

Z. = -0.029 

o 

e 

M 5 = -6.04 

c 

2 

Oj = 26.84 


! I 

: I 

j" The corresponding values of m.s. normal acceleration were: 

- T X = 1 , J = -008668 

; ’ X s 10 , J = .008973 

? .jy ■ . ■ . ■ * " 

• The percent degradation is 

X = 1 AJ = 8.55 % 

X = 10 AJ = 8.58 % 

'• J 

Thus the percent degradation is relatively insensitive to control effort 
' ' constraint. 

The degradation due to change in each individual parameter was 
calculated using the linear approximation theory, and the results are 
shown in Table 1.1; only the case X = 1 is considered. Note that the 
largest change in parameter occurs for M but the performance degradation 
V due to this change is relatively small being only 0.2%. Of course, as is 

only to be expected, the largest performance degradation was due to un- 
certainty in the control variable . It is interesting that the various 

e 

parameter changes cancel out; in particular the particular change in 

- actually causes a performance improvement ! The net degradation in 

6 

performance of the linear approximation basis yields 10.96% which should 
be compared to the actual value of 8.55% The linear approximation thus 
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tends to slightly exaggerate the degradation, at least for this particular 
case. 

It would appear that uncertainties in the control derivatives result 
in greatest degradation while the uncertainties in the stability derivatives 
play a much less significant role. Fortunately in the overall parameter 
estimation technique the control derivatives are usually the most 
accurately determinable. 


-] 9 - 
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TABLE 1.1 

Performance Degradation (Jet Star Data) 

Due to Parameter Uncertainty 

M.S. Normal Acceleration: No Control: 0.01185 

X = 1 



Nominal Parameters 

Dither Input Estimates 

APerf . 

Performance Degradatio: 


(Maneuver A) 

(Maneuver X) 


(AJ/J) 

z 

-1.45 

-1.41 

0.04 

-5.0% 

a 

i 




M 

-9.79 

-9.42 

0.37 

-0.08% 

a 





M q 

-1.43 

-0.53 

0.90 

-0.24% 

Z 6 

e 

-0.098 

-0.029 

0.069 

-7.0% 


-8.15 

-6.04 

2.11 

+1.43% 

'V 




e 





2 

a 

25.2 

26.84 

1.64 

1 

o 

o 

'■o 

on° 


Total -10.9% 





1.3 HANDLING QUALITY AND GUST ALLEVIATION. 


In this section we develop a new criterion for handling-quality im- 
provement and gust alleviation. Following usual practice, we interpret 
handling-quality as the ability of the aircraft to respond closely to 
the pilot input (signal). For this purpose we use a stochastic model for 
the pilot input, and a quadratic error criterion. We asume a Dryden or 
similar rational spectrum model for the gust, and a quadratic minimisation 
criterion as well for gust alleviation. Finally we add a soft constraint on 
the control effort. 

Mathematically, the problem can then be stated as follows. We begin 
with the linearised rigid body perturbation equations of motion of the 
aircraft: 

x(t) = x(t) + u(t) + F^ n(t) 

y(t) = x(t) + G^ n(t) 

where x(*) is the state, enhanced to include the gust generation model, 
y(*) is the senso r output , u(*) iS 'the control to be optimised, and n(t) 
is white Gaussian with unit spectral density. 


to reflect the fact that gust and sensor noise processes are independent, and 

Gf Gf > 0 (Nonsingular) . 



Let v(t) denote the pilot input . Then we can formulate the 
performance criterion as: 

T ip T 

Pr= |e f ||Lx(t) - v(t)|| 2 dt + If EC| (Qx(t) 1 1 2 ) dt + ~f E( | |u(t) 1 1 2 ) dt 
-'o 1J o 

Here L and Q are appropriately clios en rectangular matrices, and A is a positive 
constant reflecting the degree of control effort constraint desired. 


To proceed further we need to specify the model to be assumed for 
the representation of the pilot input. Ihus let: 


v(t) = C 2 z(t) + G 2 N(t) 
z(t) = Aj z(t) + F 2 M(t) 

where M(t) is white Gaussian with unit spectral density, with 


F 2 g“ 


G 2 G 2 > 


0 

0 (Nonsingular) 


In what follows we shall specialise to 


G„ G* = cl , I = Identity matrix. 

2 2 ’ * 

where 

e > 0 

as being general enough. If e>0, this is tantamount to saying that we alia-: for 
seme error in the sensor ; we shall consider the case where e is zero, separately 
below. 


The key step in our theory is to recast the performance .. 
criterion, making essential use of the pilot inrut rcdel, into the •'stochastic 
regulator" form., from which the optimal feedback control can be readily deduce.'.. 



It is convenient to specify the dimensions of the various 
vectors and matrices we shall use to clarify the notation . Thus let: 


x ~ n x 1 

z ~ m x 1 

~ n x j 
C, ~ r x n 

~ p x m 
L ~ p x n 

Q ~ q x n 

Let X denote the (n+n) x 1 matrix: 



Then 


| |L x - C ? z| j 2 + 1 |Q x 1 1 2 


where 


L"L + Q\> 


Tv. R X X 

A 

-L" C, 

c" c„ 


Note that R is (m+n) x (n+n) and is of course ncn-nerotivv definite. 

Using, XCt) as the conreund state space, and R as do 
above, the perfoiTimce criterion Bin can no..' N* rewritten as: 





Pr. = (1/T) f El[R X(t), X(t)] + A[u(t) ,u(t) ] ) dt 

J 0 

where C , 3 denotes ii'ner product, and where X(t) satisfies the ’state' equation 

x(t) = Axct) + BuCt) + F n<t> 


where 


A 


B 


A x 0 
0 


B, 


(n+nOxj 
0 


n Ct) = 


n(t) 

Il(t) 


The control u(t) is no.-? to be based on the "observation" vector: 

yCt) 

_ v(t) 

where of course: 


y(t) 

v(t) 


c, o 

x(t) 

G. 0 

1 

0 c. 

+ 

c(t) 

1 

0 . G 

/ 


C 


n(t) 

!i(t) 



where 


The optimal feedback control u Q (t) can then be expressed as: 


u (t) = (-1/A) B* ( pj; 1 xCt) + P 12 l(t) ) 

o 1 c c 


Finally 


and 


x(t) 

= A 1 x(t) + 

P^ C* ( y(t) - 

C, x(t)) + B u (t) 
1 o 

i(t) 

= ^ z(t) + 

2 * 

P* C ( 
f f u 2 t. 

v(t) - 

C 2 z(t>) c -1 

p*(t) 

= A 1 pJ(t) ■ 

^ P^Ct) 

t - fJ 1 ' 10 ! C 1 

F^(t) 

= A 2 P 2 (t) • 

^ P 2 (t) 

ft ft ? _ ft 

A + r r _ rr(t) C C 

*7 *2 *2 *r X; L 2 2 



> 




P 





c 

/o 

22 

P Z/ 

c 


(t) + 

A* p (t) + p 

' C ( 

,(t) A - 

P c (t) B B* P r (t) + p = o 


-1 


As in most applications, we may use the ' stead;,’ state' versions of the sclut 
to tlie Pdccatti equations. 

Either the handling quality criterion or the [rust 

alleviation criterion ( or both) can be expressed in the f err: 

,T* 

tJ ll:: x(t)ll l dt 


which can be calculated as: 



• . < 


Tr. M* M P f + Tr. M* M J 
where J is the solution of the linear equation: 

0 = (A - B B* P c |A) J + J (A - B B“ P c |A) + P f C*C P f 


where 





Example 

As a specific example we shall take the one-dimensional pilot model: 


z(t) = A^ z(t) + F^N(t) 

where 

Ay = -(Vl.ll) 

F 2 = Identity 

we shall also take the observation to be noise free 


[e=0]. Ln this case 



[and more generally where the full state is observed noise-free, or C 2 is 
non singular and e=0] 

z(t) = z(t) 


and the combined system of equations for x(t) and z(t) is given by: 


.1 * 


x(t) = A 1 x(t) + C 1 (y(t) - C.jX(t)) + B^u^Ct) 


z(t) = A 2 z(t) + F 2 N(t) 


We shall also specialize to the case where we are not interested in gust 
alleviation but only in handling-quality so that 

Q = 0 

We shall further take X = 1 for simplicity. 

11 12 

The main thing is the calculation of P and P . For this we must 
° c c 

solve the Steady State Riccatti equation for P c which in the present case 

is (dropping the subscript c) : 

ft n 11 n ft n ft 

A 1 P + P A x - P • l B 1 B 1 P x± + L L = 0 

21 * 91 91 ft ft 

P A 1 + A 2 P " P Vl P H " C 2 L = 0 

ft 22 22 21 * 12 * 

a 2 p“ + p a 2 - p b 1 b 1 P 1 ' + c 2 c 2 = 0 

°2 01 11 

Because A 2 is 1 x 1 we can express P" and F" in terms of F . Thus we ra 


21 1 
P X (A 1 + A 2 


t_b B P ) 

- i i *ir 


ft 

= C 2 L 


-26a- 



or 


ft 23 * * ft _i 

P 12 = P = C 2 L( *l + A 2 - B 1 B 1 P ll> 


22 * _1 / 21 ft I? ft 

P = (A 2 +A 2 > P Vl P ’ -C 2 C 2 


Because e=0, J is now given by: 


* ft ft 

(A - BB P )J + J(A-BB F ) + 

c c 


r 1 * 1 

P f C l C l P f ° 


f 2 f 1 


= 0 


And finally the mininal value of 


lim 

T-xo 


T f 


E| |L x(t) - v(t) j | Z dt 


Tr H H J + Tr. M li' 


r P* 0 
0 0 


where 


M M 


* * 


L L 

-l c 2 

* 

- C 2 L 

c"c 


Decomposing J as: 


J = 


J 11 J l? 


‘"'21 
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we have: 


*2 J 22 + J 22 A 2 + F Z F *- = 0 


(Ai^s/'P 11 ) J 12 - RftV 2 J 22 * J 12 a/ = 0 


(Ai-B^V 1 , J u - ^\V \ 2 


* J 11 C V E L B 1 V1) ' - J 2 2 p21 Vl‘ 


1 * 1 

+ c n c, ?i = o 

X 1 -L X 


Taking advantage of the fact that in our example A 2 is one-dirrensional , 
we have 


J 22 = -CA 2 +A^) -1 F 2 F* 


etc. , and thus expressing J-q> J-^ 2 ^ terms ^ 22 ‘ 

Finally 


Tr. M M*J = Tr. (L LJ n - L *C J 21 - C 2 *L J u + C" C ? J„) 


Taking C 2 = 4 and the longitudinal Jet-Star equation at the nonunal values 
(p. 57) we have: 
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\ . - 

'Hie control gain matrix P is calculated to be: 

c 

1.32 -1.63 -.224 

-1.63 2.02 0.277 

-0.224 0.277 .033 

0 0 0 

.0811 -0.1 -.0138 

P* 2 = (10 -7 ) [4.935, -6.11, -.8393, 0, -1.371 x 10 6 ] 

22 

P = 4.710 

c 

-b’ ! p 11 = 10~ 7 [-1.622, 2.008, .2759, 0, (-20) (2.07) 10 5 ] 

1 c 

Hence we see that with good approximation 

u o (t) = (-20) (2.07) 10' 2 x 5 + (20) (1.37)10' 1 z 

= (20) (10 -2 ) ( 2 . 07 £ 5 - 13.7 £) 

where = 6 f (t) , the actuator output. As expected, the controller virtually 
ignores all other states. Note that the control gains arc independent of 
the gust intensity. 
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The corresponding value of the steady-state mean quadratic deviation 

= Tr. MM”j + Tr MM*P f 

= 2.58 +0 

= 2.58 

The mati’ix J is 


2.13 

5.37 

3.01 

o 

-1.62 

.0 

5.37 

4.59 x 10 4 

4.59 

0 

-5.07 

-2.03 

3.01 

4.59 

9.40 

.0 

-3.64 

-1.83 

0 

0 

0 

0 

0 

0 

1.62 

-4.07 

-3.64 

0 

2.02 

1.04 

0.08 

-2.03 

-1.83 

0 

1.04 

0.0556 



IDENTIFICATION OF DERIVATIVES FROM FLIGHT DATA: THEORY, ALGORITHMS 
AND SIMULATION RESULTS. ... 
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In this section we report on Task B: development of an algorithm for 
extracting stability and control derivatives from flight data taking gust 
into account. Section 2.1 deals with the essential background theory. 

The algorithm itself is developed in Section 2.2. Section 2.3 details of 
verification of the algorithm on simulated flight data. 

2.1 THEORY 

The theory will be presented in a degree of generality that car. handle 
the linearised model of any rigid-body aircraft motion, longitudinal-and- 
lateral, and taking gust loading into account. 

We use a ’ cont inous -t ima ’ model rather than a 'discrete-time model 
(as in most of the engineering literature) primarily because the data 
sampling rate (50 samples/sec) is much higher that the Nyquist rate cor- 
responding to the noise bandwidth [5 Hertz at most}, so that the usual 
assumption of independence of noise samples on the output is not valid. 

On the other hand, time-continuous model 'integrals' are well approximated 
by 'sums'. A fuller discussion of the trade-offs in the choice between 
the two models my be found in Appendix I. 

It is convenient to rewrite the motion equations and sensor data in 
state-space form. Thus we have: 

x(t) = A x(t) + E u(t) + F n(t) 

for the 'state 1 equation, where the matrix A contains the stability 
derivatives, and B contains control derivatives, u(t) being the input 



(6 6 6 etc.) and n(t) is white Gaussian noise with unit spectral 

e a 5 r 

density matrix. The term F n(t) characterises the gust, as we shall see 
below. We denote by the column vector Y(t), the various sensor outputs 

(observed data) and can be expressed as 

Y (t) = CX(t) + D u(t) + G n(t). 

The 'state' noise F n(t) and "observation noise" G n(t) are independent: 

F G* = 0 

and 

GG = diag. [d^3 

with all positive - corresponding to the fact that we do not assume 
that any sensor is noise-free. 

We denote by 0 the (vector of) various unknown parameters in all the 
matrices A, B, C, D, F, we estimate 'O' by minimising with respect to 0, 
the unknown initial state and possibly unknown noise matrix diag. 

| d^ | , as well as bias : 

q (0;T) 

= Tr. Log GG" + 1 f C(Y(t) - Cw(t) - D u(t)), 

T J 0 

(GG*)” 1 (Y(t) - Cnft) - Du(t))1 d 

+ Tr. ( GG *)" 1 C F C* (2.1) 


_~Q_ 



where x(t) is the Kalman state estimate defined by: 


(t) = CA - P C*(GG*) _1 C) x Ct) 


+ B uCt) 

+ PC* (GG*)** 1 Y(t) 


X(0) = Xr 


and P is determined by: 


( 2 . 2 ) 


* f: ft ft -1 

o = AP + PA + FF - PC (GG) CP 


(2.3) 


Let us next calculate the C - R bound matrix, assuming that (GG ) 
is known. Let 

y (t) = Cx(t) + D u(t) 


R (0 ;T) = 


}e f [ (GG*) 1 i V (t), S Y (t) ] dt> 

* J ° V I 


(2.M ) 


where |ou| denote the components of 0. Then letting 0^ denote the true 
value of 0, the C. R. bound matrix is 


R (0 Q } T) 1 
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where 0 is the true -value of 0 (no model error), 
o 

we approximate this by 

(T R (0 T ; T) r 1 

1 

dt (2.5) 

* 

Note that the approximation depends on the observed data directly (in 

A 

addition to the dependence via Q^, the estimate). 

2.2 ALGORITHMS 

In this section we develop the algorithms for parameter estimation in 
the generality of the formulation of Section 2.1. 

In minimising (2.1) we first use nominal values for the stability 
derivatives, and take G G as the identity matrix, unless there is a 
priori values of | available from previous flight test data in clear 
air. In the first iteration step we vary only the Control derivatives 
(in B and D) as well as the biases. Note that in this process A, C, 7, and 
P remain fixed. A first estimate for (G G ) is then obtained using: 


where 0 is our estimate for G and 
T 


R ( 0; T ) 


'{* / 


[ (G G*)” 1 3_ Y (t), 3Y (t) ] 

0 3 3 

a. a. 

i 1 



1 / (Y (t) - Y (t) )? dt 

T J 0 


where ( )„. denotes i component. 

using this value of (G G ) we now iterate with respect to all the parameter 
0, c, and the biases, using the following formula: 




0 r - R (O n ;T) _1 G CG n ;T) 


where G (0;T) is the vector with components 


/ T f /s 

[(GG V 1 3 Y Ct), YCt) - YCt) J 

0 T 

T 


dt 


+ Tr. (G G*‘) -1 3 CP c' 

V 

l 


and R (0;T) as defined already in (2.5) 

A few words on the details of the computer calculation of relevant 
quantities may be in order. The main step is the calculation of the partial 
derivaties of dynamic variables: 


| xCt) 
a. 

l 


For this, we take the first order- variation on the Kalman equation: 


H\ iM ) 


= (A -PC 03 >:(t) + 3 (A - ? C" C) S(t) 


+ / 3 BVu(t) 
3' 


+ (l PC 1 

' nr . ' 


(G G ) _1 YU) 



time continous linear 




where A is the sampling interval CQ.02 sec. in our case). 

The Riccati matrix P is computed by the following algorithm: 
[we assume A-Stable , A - F controllability and C - A observability} 
compute A from: 

A" A + A A + C* C = 0 

then solve (the linear equation) 

0 = (A - P C" C) P + P.. (A - PC* C) " 
n n+x n+X n 

+ F F" + P c“ C P 
n n 

with 




2.3 VERIFICATION OF ALGORITHMS : SIMULATION RESULTS. 

The algorithms developed in the previous section were verified by 
simulating the longitudinal short -period motion dynamics. The equation 
used in State Space form were : 

I xCt) = Ax(t) + Bu(t) + Fn o (t) 
v(t) - Cx(t) + Du(t) + GnCt) 

where 


— — 




a 


0 


s 


m 






r *1 

0 


0 


5 



m 


e 

• 





0 

v = 

a 

u = 




n 


_ 1 



m 



a 




g_ 


a 



- m - 







— 




— 

z 

a 

. 0 

1 

z 

a 

. 

Z 6 

e 

z o 

0 

0 

1 

0 

E = 

■ 

0 

0 

K a 

0 

0 

K 

a 


V 

“5 

e 

K ° 1 

0 

0 

0 

-v/1000 _ 


- 0 

0 - 1 





and n(t) and n g (t) are scalar and vector valued independent Gaussian v.’hite 
noise processes with zero mean and unit spectral densities, respectively. 

We assure that G is known, and estivate K a ,!i-, Z . , end 

a a o Q c- e 

2 

o from sensor measurements v(t) and tacv.n control inrut 6 (t). 
vj e 

g 
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Input valves for the simulation are listed in Table 2.1. The control 
input <$ e is taken to be a square wave of amplitude 0.02 rad. and of frequency 
0.4 hertz. 

The Power Spectral Density (PSD) based on 512 points of simulated 
measurement noises are shown in Figure 2.1. The horizontal line is the 
desired spectral density assuming the noise process in white. The noise 
processes shown in the figure are seen to be reasonable approximations of 
white noise. 

Figure 2.2 shows the PSD of a , and its theoretical asymptote. Since 

6 

the turbulence is direct proportional to a^» its PSD will also have the 
same shape. 512 points were used in the computation. 

The rate of convergence of the Riccati matrix, P, to its steady value 
is examined. The diagonal elements are plotted in Figures 2.3 and 2.4. 

Note that P converges to its steady state value after approximately 12 
iterations . 

Figure 2.5 is a comparison of the simulated time -history and the 
estimated time-history obtained from the identification algorithms. The 
solid lines are' simulated data and the broken lines are the estimated data. 
The fit is observed to be excellent. 

The residuals are plotted in Figure 2.6 and their FSD's are plotted in 
Figure 2.7. The residual PSD's are essentially flat over all frequencies. 

The convergence of the estimated parameters based cn 5.12 seconds of 
data shown in Tabic 2.2. All parameters converged in sic: iterations. Th- 
final estimates of all parameters with the exception of the last o.ne arc 
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2 

very close to the true parameter values. Since the last parameter a is 
the turbulence power, it should measured in decibels. Here the difference 
is less than 2 db, 

Finally, the identification algorithm is applied to other simulated 
data records. Each record has the same control input but the state and 
measurement noises are independent from those of the record. The mean and 
variance of estimated parameters are then computed. Tables 2.3 thru 2.G 
mean and variance with the true parameter value and the Cramer-Rao bound. 
The mean values are very close to the true parameter value. The variances 
are reasonable close to the Cramer-Rao bound, for the small number of 
samples . 
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Table 2.1 


Data for Simulation 


- 1.65 rad/sec 

- 54.0 rad/sec 

- 1.65 rad/sec 

- 0.45 rad/sec 

- 52.5 rad/sec 2 
25.0 ft 2 /sec 2 

Z o= 0 

K = 0 
o 

V = 1670. ft/sec 

l = 10 ft 
z 

l = 32 ft 
a 

K = 1.7 
a 

A =0.01 sec 
g^ = 0.0005 rad/sec 
g 2 = 0.0001 rad/sec 
g 3 = 0.01 g 

g 4 = 0.00005 rad 




r 


Table 2.2 Convergence of Estinvated Parameters Using 5.12 Seconds of Data 


ITERATION 

2 

a 

M a 

M. 

6 

Z 6 

e 

M 6 

e 

a 2 

w 

g 

Initial start- 

- 2.40 

- 39.0 

- 2.40 

- .675 

- 36.0 

2.5 

up value 







1 

- 1.764 

- 53.2 

- 1.887 

- .675 

- 51.09 

3. 85 

2 

- 1.646 

- 53.9 

- 1.655 

- .447 

- 52.42 

6.95 

3 

- 1.660 

- 54.3 

- 1.692 

- .464 

- 52.88 

13.06 

4 

- 1.650 

- 54.0 

- 1.665 

- .446 

- 52.56 

13.47 

5 

-1.646 

- 53.8 

- 1.647 

- .445 

- 52.31 

19.98 ' 

6 

- 1.641 

- 53.9 

- 1.651 

- .446 

- 52.50 

20.05 

True value 

- 1.65 

- 54.0 

- 1.65 

- .45 

- 52.5 

25.0 
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Table 2.3. Statistics of Estimated Parameters 
No. of independent runs = 10 

No. of data point for each run = 256 (2.56 seconds) 

SUfJLATICX: TEST 


PARAMETER 

TRUE VALUE 

C-R BOUND 

MEAN 

VARIANCE 

Z a 

- 1.65 


- 1.650 




.00014 


.00022 

N a 

- 54.0 


- 54.07 




.057 


. 066 

M. 

- 1.65 


- 1.636 . . 


6 







.00079 


.00072 

Z 5 

- 0.45 


- .4525 


e 







.00016 


.00018 

m 6 

- 52.5 


- 52.55 


e 







.070 


.057 

2 

a w 

• 25.0 


18.76 


g 







22.9 


4.13 
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/ 


I 

I 

* i 

j 

Table 2.4. Statistics of Estimated Parameters 

" No. of independent runs =18 

No. of data points for each run = 512 (5.12 seconds) 





simulation 

TEST 

PARAMETERS 

TRIE VALUE 

C-R BOULT) 

MEAN 

VARIANCE 

Z 

a 

- 1.65 

.000045 

- 1.650 

.000042 

M 

a 

- 54.0 

.019 

- 54.01 

.026 

K. 

6 

- 1.65 

.00027 

- 1.6-.1 

.00031 

Z 6 

e 

- 0.45 

.000057 

- 0.4504 

.000035 

H 6 

e 

- 52.5 

.025 

- 52.51 

.025 

o 2 

w 

g 

25.0 

6.7 

17. *.3 

2.5 
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I 

I 

Table 2.5. Statistics of Estimated Parameters 

! 

f No. of independent runs = 10 

No. of data points for each run = 1024 (10.24 seconds) 

r 

SIMULATION TEST 


PARAMETER 

TRUE VALUE 

C-R BOH© 

MEAN 

VARIANCE 

Z 

a 

- 1.65 


- 1.651 




.000021 


.000015 

M 

a 

- 54.0 


- 54.03 




.0082 


.0093 

M. 

- 1.65 


- 1.650... 


6 



* 




.00012 


.00014 

Z 6 

- 0.45 


- .4502 


e 







.000026 


.030029 

h 

- 52.5 


- 52.55 


e 







.011 


.011 


25.0 


18.48 


g 







3.28 


1.66 
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Table 2.6. Statistics of Estimated Parameters 


No. of independent runs = 8 

No. of data point for each i\m = 2048 (20.48 seconds) 

SIMULATION TEST 


PARAMETER 

TRUE VALUE 

C-R BOUND 

MEAN 

VARIANCE 

Z a 

- 1.65 

.0000082 

- 1.650 

.00001 

M a 

- 54.0 

.0032 

- 54.03 

.0029 

M. 

e 

- 1.65 

.000046 

- 1.650 

.000072 

z 6 

- .045 


- 0.4499 


e 


.000011 


.000015 

M. 

- 52.5 


- 52.54 ' 


e 


.0044 


.0023 

1 

a 

V 

25.0 


18.42 


B 


1.29 


2.94 
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3. FLIGHT TEST DATA: 


LOCKHEED JET STAR 



3. 


FLIGHT TEST DATA: LOCKHEED JET STAR 


In this section we evaluate the performance of the identification 
algorithms on flight test data taken on a Lockheed Jet Star flying in 
turbulence (courtesy of K. Iliff, DFRC). 

3.1. STATE SPACE MODEL DYNAMICS 

The Flight data contained only the longitudinal short period mode. 

The turbulence was modelled by a Dryden spectrum, in which both the 
intensity and break frequency were allowed to be parameters to be estimated. 

The sensor data vector v(-) consists of: pitch rate, pitch, nomal 

acceleration, and angle of attack. The state equations can be written: 

x(t) = A x(t) + B u(t) + F n (t) 

s 

v(t) = C x(t) + D u(t) + G n(t) 

where the state (enhanced to include the gust cc^pcrent) is the column 
vector : 


angle of attack 
pitch angle 0 
pitcn rate 


0 



A 


2a 
0 

M a 
0 

Z 6 

e 

B = 0 

]I 5 

e 

0 
6 

u = e l 

The rest of the iratrices are as in section 2 , except for matrix G of 
error densities: 

G = Diag. ( C.0005, 0.0001, 0.01, 0.00005) 

In the control derivative matrix B, it is optional whether Z. and 

6 e 

M. are allowed to be independent parameters or not. It would appear 
c - - 

e 

that making them dependent (one a multiple of the other) yields better 
fit. 

The time history data was subdivided into five segments of sup-maneuver 
marke A, B, C, D, and E- 

Each segment was processed independent 3 v of the others. The 
segments had to be combined to study the effect of changing the tor science 


0 

0 

0 

0 


1 
1 

Ma !1 


a 

-V /1000 


0 

0 

0 
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spectrum break frequency. Table 3.1 lists the nominal Lockheed Jet Star 
data. 


3.2. EXPLORATORY COMPUTER RUNS 

Prior to actual processing, many exploratory computer runs were rede 
to determine the relative role of the various parameters auxiliary to the 
min derivatives. 

3.2.1 ITERATIONS REQUIRED TO MINIMISE COST FUNCTIONAL: TYPICAL CASE. 

In addition to the parameters already listed, additional bias 

parameters as well as initial states were included. Figure 3.1 shows 

the behavior of the cost functional as a function of the number of 

iterations. As my be expected it drops steeply in the beginning and 

flattens out after five or six iterations. The maneuver chosen is 

maneuver A; the wind tunnel values of the derivatives were chosen for the 

starting values, and the initial states all set to zero. Table 3.2 shews 

the actual values of the cost functional for each iteration. The initial 

sharp drop occurs as a result of adjusting the 'lineary dependent' 

parameters: the biases and the control derivatives. A plot of the 

behavior of the derivatives in shown in Figure 3.2. The gust scale-lengtr. 

was taken to be L = 1000, the intensity being allowed to be unknown. 

W £ 

The starting value for 



was taken to be zero. 
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3.2.2. EFFECT OF CHANGING STARTING DERIVATIVE VALUES. 


It is of interest to determine the effect of choice of starting 
values of the derivatives on convergence of iteration and the final values 
of the derivatives estimated. This is illustrated using maneuver A. 

Table 3 shows the values of the derivatives for various choices of the 
starting values, as well as the corresponding final (minimal) values of 
the cost functional. The starting values of the biases and states were 
set to be zero in all cases. Note that the final values of derivatives 
are relatively insensitive to the starting values. The number of iterations 
to attain the minimising cost functional value were always four or five, 
and not sensitive enough for inclusion in the table. 

3.2.3. EFFECT OF DELAYING ANGLE OF ATTACK (MEASURED). 

Because of the physical displacement cf the vane measuring angle 
of attack, an attempt was rede to see whether delaying the angle of attack 
time history with respect to the rest of the data would have any effect cr 
derivative extraction. This is shown (Run 6) in Table 3.3. It is conduces 
that no significant improvement results from doing tills . In the rest cf 
the work reported the angle of attack data was not delayed. 

3.2.4. EFFECT OF CHANGING BREAK FEZCUE.CY IN DRYDEK SPECIE JK. 

To determine the effect of charging the (3-db) break frequency in ft 
gust (Drycer) spectrur . maneuver A was run (with starting values set at t'n 
wind-tunnel values), clanging tne break frequency to values an octave ab-.t 
and an octave below. As can be seen trere is virtual Iv no effect on 
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T A B L E 3.1 

LOCKHEED JET-STAR DATA (NO! DUAL) 

V = 600 ft/sec 

Weight = 3500 lbs. 

Sampling interval: .02 sec. 

= 0 

l = 26 ft 
a 

K = 1.6 
a 

Z = -1.5 rad/sec 

a 

I! = -8.82 rad/sec 2 

a 

Mg = -1.4 rad/ sec 

Zg = -0.11 rad/sec 

G 2 
= -10.9 rad /sec 

e 



maneuver A, as is to be expected because the time history duration is rot 
large enough. To remedy this, maneuvers A, B, C, and D were combined, and 
same changes made in the break frequency. Figures 3.3 and 3.4 show the 
estimated gust spectrum. The effect even here does not appear to be 
significant, because the time-history is still not long enough. 

3.2.5. EFFECT OF FAKING Z_ DEPENDENT ON M. 

6 0 
e e 

As mentioned earlier, the control derivatives and M. ray be 

°e °e 

allowed to be independently determined, or one allowed to be a multiple of 
the other, and the coefficient of multiplication estimated. Table 3.3, 
run 9 shows the effect of making them dependent. As compared to the 
standard run, run 1, in the same table, there does not appear to be any 
significant change, except perhaps in the value of itself. 


3.3 MAIN RESULTS 

In this section we present the rain results of derivative extraction. 

The results on each maneuver are grouped arid identified by the segment of 

time-history: A, B, C, D, E, CD, and ABCD. In all the cases, the starting 

values of the derivatives were the wind-tunnel values, with the initial 

states all set to be zero, ans L taken to be 1000. 

w 

£ 

For each maneuver a set of five figures, each clearly identified by 
the segment used are given the tire history plot superimposed or t’->c 
estimated tima history for visual cor - arisen , the tirr history plot cf r 


errors’ 



Aq 

A0 


Aa 

n 

Aa . 

and the corresponding P.S.D's and 
estimated states: 


0 

A 

q 


a > 
g ’ 


Note that a yields (an estimate of) the actual gust incident waveform. 

& 

The p.s.d. derived from this, the estimated gust p.s.d., is also shorn. 

The figures have been numbered consecutively, from Tig. 3.5 thru fig. 3.39. 
The parameter values for all these maneuvers and the C.P. bounds (in 
parenthesis) are tabulates in Tai le 3.4. Finally figure 3.40 shows the 
means and variances over these maneuvers, along with the wind-tunnel 


values. 






TABLE: 3.2 


Iteration 

Cost Functional 

0 

324558.27 

1 

2303.46 

2 

-918.270 

3 

-1417.5673 

4 

-1493.9562 

5 

-1497.4350 

6 

-1497.5041 

7 

-1497.5119 

8 

-1497.5084 

9 

-1497.5103 

10 

-1497.509*1 

11 

-1497.5099 

12 

-1497.5097 

13 

-1497.5099 

14 

-1^97.5093 

15 

-1497.5099 



TABLE : 3.3 


RUN ' 



M 0 


% 

B 

Z 0 

M 0 

a s (0> 

q(0) 

a E (0) 

Cost Funct. 

1 

-1.45 

-9.79 

-1.43 

-.098 

-8.15 

25.20 

.011 

-.276 

-.015 

-.0022 

.015 

-1497.1992 

2 

-1.45 

-9.79 

-1.43 

-.098 

-8.15 

25.18 

.011 

-.276 

-.015 

-.0022 

.016 

-1497.2302 

3 

-1.45 

-9.79 

-1.43 

-.098 

-8.15 

25.17 

.011 

-.276 

-.015 

-.0022 

.016 

-1497.2401 

4 

-1. a 5 

-9.79 

-1.43 

-.098 

-8.15 

25.19 

.011 

-.276 

-.015 

-.0022 

.016 

-1497.5099 

5 

-1.45 

-9.76 

-1.44 

-.096 

-8.17 

25.78 

.009 

-.277 

0.0 

0.0 

0.0 

-1497.2622 

fi 

-1.46 

-9.86 

-1.57 

-.066 

-8.15 

26.95 

.012 

-.276 

-.014 

-.0006 

.013 

-1498.8137 

7 

-1.46 

-9.89 

-1.43 

-.110 

-8.24 

11.95 

.010 

-.279 

-.010 

-.0021 

.011 

-1497.4225 

n 

0 

-1.45 

-9.79 

-1.43 

-.098 

-8.16 

50.13 

.012 

-.276 

-.027 

-.0022 

.028 

-1497.2538 

0 

-1.44 

-9.73 

-1.44 

-.075 

-8.04 

25.90 

.012 

-.272 

-.015 

-.0021 

.016 

-1496.9238 


Start at Wind Tunnel Value. 

Start at -1.3, -10.84, -1.83, -.001, -7.16,10.0, 0, 0, 0, 0, 0, 0, 0 

Start at -1.0, -15.0, -2.0, 0., -5.0, 2.0, 0, 0, 0, 0, 0, 0, 0 

Start at -2.0, -5.0, -.5, -.3, -15.0, 2.0, 0, 0, 0, 0, 0, 0, 0 

Start at Wind Tunnel Value, x(0) = 0. 

Start at Wind Tunnel Value, a measurement delayed .04 second (2 data points) 
Start at Wind Tunnel Value, =500 
Sturt at Wand Tunnel Value, L v , = 2000 

Start at Wind Tunnel Value, L Set Z. =F M. , F =0.00932 

’ c <5 c 

e e 
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TABLE 3.4 Composite Derivative Values 



Z a 

H a 

M 0 

Z 8 

e 

M 6 

e 

g 

A 

-1.45 

(.0472) 

-9.79 

(.349) 

-1.43 

(.1297) 

-.098 

(.0355) 

-8.15 

(.329) 

25.20 

( 5.3) 

2 

-l. u 6 

(.0505) 

-9.79 

(.380) 

-1.48 

(.1233) 

-.0077 

(.0437) 

-7.35 

(.404) 

63.51 

( 12.6) 

c 

4 

r ~* 

t 

(.0429) 

-9.66 

(.321) 

-1.86 

(.1039) 

-.150 

(.0335) 

-9.00 

(.294) 

75.00 

( 13.2) 

1** 

-:.“j 

(.0449) 

-9.76 

(.330) 

-1.18 

(.1006) 

-.011 

(.0331) 

- 6.68 

(.302) 

65.04 

( 13.5) 

- 

-1.31 

(.0315) 

-3.03 

(.220) 

-1.83 

(.0788) 

-.044 

(.0310) 

-8.03 

(.244) 

305.98 

(153.7) 


-1.4U 

(.0337) 

-9.89 

(.255) 

-1.31 

(.0826) 

-.031 

(.0275) 

-7.80 

(.251) 

45.22 

( 6.4) 


-l. ; 7 

(.0312) 

-10.01 

(.239) 

-1.62 

(.0747) 

-.078 

(.0230) 

-7.76 

(.208) 

71.76 

( 9.3) 

• 

-1.45 

(.0225) 

-9.94 

(.174) 

-1.43 

(.0537) 

-.062 

(.0175) 

-7.65 

(.159) 

57.70 

( 5.3) 

- 

-1 . 

-7./0 

-1.31 

-.123 

-9.75 
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FLIGHT TEST DATA: 


EVALUATION OF DITKER-L’ JFUT 


4. FLIGHT TEST DATA: EVALUATION OF DITrZR-IIJPUT. 


One of the unresolves questions in parameter identification is the need 
for a conscious manual input by the pilot. In all the identification 
cases studies and reported, time histories were chosen which contained an 
impulse-like pilot input (cases A, B, C, D, E). In order to study the 
possibility of using a dither-like input instead, we took portions of the 
data where there was no sustained pilot input. Instead of the segments 
A, B, C, D, E as before, we used the segments marked "X' 1 and "Y" on 
Figure 4.0 which shows the total Jet Star tine-history available. The 
section X and Y are well beyond the pilot impulses, the pilot input 
excursion being now almost "dither-like". 

The parameter extraction program was applied to the segments X and 
Y, and the results obtained are shown in a series of figures. Figure 4.1 
shows the time -histories of the sensor data as well as the (dither) inputs 
corresponding to the segment X. The estimated values for G, tu, n , and a 
are shown in broken lines in the sa^e figure for comparison. Since the 
scale chosen may tend to obscure the discrepancy. Figure 4.2 shows the 
actual error time-histories, now to a differenct scale. Figure 4.3 show? 
the (Kalman) estimates of the states. The p.s.c. of tne estimated 
turbulence is given in Figure 4.4, and tre p.s.d.'s of the fit errors 
are shown in Figure 4.5. The corresponding quantities for maneuver i are 
given in Figures 4.6 thru 4.10. 

Hie estimated parameter values for bet), dither segments, as .ell at 
the contiguous pilot mar -invert / and L, are t~o; n in Table 4.1, wheie • 


J » 



have also shown the wind-tunnel values for purposes of comparison. The 

large error in ML is expected, being a control derivative. The discrepancy 
‘’e 

in M^, a stability derivative, is a surprise. The values of and 

obtained are close to those obtained in maneuvers A and E, as expected. 
The corresponding degradation in control performance has been calculated 
in section 1. 
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Parameter Segment X 


Maneuver A 


Segment Y Maneuver E 


z 

a 

- 1.41 

-1.45 

-1.38 

-1.38 

M 

-9.42 

-9.79 

-7.34 

-8.03 

a 





M e 

-0.53 

-1.43 

-1.31 

-1.83 

H 6 

e 

-6.04 

-8.15 

-10.49 

-8.03 

9 

0 

26.84 

25.20 

1177.5 

305.98 


TABLE 4.1 

Estimated Parameters: Dither Input and Gust 
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APPENDIX I 


In this Appendix we examine a question basic to all digital computer 
processing of flight test data: whether the basic aircraft dynamics should 

be modelled as a continous-time process or a discrete-time process, tie 
show that even though all digital computers sample the data before processing, 
the choice between the two models - continous-time vs. discrete-time is not 
clear-cut - that, in fact, for flight test data processing considered in 
this report, the continous-time model is more efficient, overall. 

It is sufficient, for this purpose, to study the linear (Kalman) 
filtering problem. In order not to clutter the main thread of ideas with 
unnecessary mathematical details we shall only consider the* case where the 
state-space dimension is one and there is only one observation. Thus let 
the system model be: 


x(t) 

= a x(t) 

+ f n (t) 
s 

y(t) 

= c x(t) 

+ g n(t) 


where n (t) is a white noise process with unit spectral dens it, - . The 
s * r 

3db signal bandwidth (two-sided) is then 

2_a 
2 - 

Vie take a = -1, typical of tno longitudinal short -per: od mode frequent . In 
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continous-time Kalman filtering theory we take n(t) to be white Gaussian. 
Let the spectral density of n(t) be unity, and let us take c = 1. Then 
the steady-state Kalman filter is given by: 


X (t) = (a - p) x (t) + p y(t) 


where p is given by 


P = “5 




2 2 
g4 + fg 


The main difficulty in applying the continuous-tire theory is that 'g' is 
essentially unknown, since ’white noise’ is only a fiction, and what is 
really meant is that the noise -bandwidth is large compared to that of the 
signal. Now, if the time-history is long enough, we ray allow 'g' to be 
a parameter to be chosen by 'minimising' (since we are taking c = 1): 


1 

T 



- x (t)j 2 dt 


We nay also estimate analytically the quantity: 


E I(x(t) - x(t) 2 J 


as a function of 'g' for assumed form of the spectral density' of the noise 
n(t). Our first model for the spectral density is tne par d-Iimitoc ve*‘s: , 
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Case 1 


p (X) = a_ -B < X < B 
S 1 2B 

= 0 


since this is not 'physically realizable', we choose next the more 
realistic : 

Case 2 


4 u f c a 2 
(2ttX) 2 + (27Tf c ) 2 

Here the 3-db break frequency is f Note that in eitner case the total 
power is a 2 . Since 

x(t) - x(t) = e(t) 


P c (X) 
S 2 


satisfies : 


e (t) = -(1 + d/£ 2 ) e (t) - f n Ct) 

* s 

+ (p/f 2 ) n(t) 


we can calculate tnc filtering error 



E Ee(t) 2 : 


by integrating the spectral density: 


for Case 1 


2 E_ 

f + ^ 


2 

? 

sF 


3 1 (A) 


4tt 2 A 2 + (1 + pc 2 /g 2 ) 2 


and for Case 2 


2 — 
f + g2 

4tt 2 A 2 + 1 + (p/g 2 ) 


4n f 0 
c 


2 2 2 2 
4tt A + 4tt f c 


The results of calculations are shown in figures A-l and A-2, 
R c versus g. Note that the optimal choice of g for Case 1 is 

2 2 
g = o_ 

2B 


while for Case 2 it is 


2 

o 



2f 

c 


as may be expected intuitively more or less. 


which plct 
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Let us now consider the situation we really want to study where the 
observation y(t) is sampled with sailing interval A. In most engineering 
treatments, we then proceed to write down the corresponding discrete-time 
version of the system model : 


x = e A x + T N (nA) 
n+1 n s 

where N g (0 is a white noise (sequence) and 


y^ = x„ + n (nA) 
n n 



E E n(nA) 2 J = a 2 

E [ N (nA) 2 ] = 1 

s 


In the discrete-time Kalman filter the assumption is made that the noise 
samples n(nA) are uncorrelated, leading to the filter equations: 


Vl = (1 ' k) iA *n + k y r 


- 12 ? 



where 


k = a 2A P + r 2 


e 2A p + r 2 + f 2 


P = e 
2 


2A 


(■ ?t V 


I 2 * nr 2 e’ 24 


° 2 ) 


„ r 2 2 -2A . 

C = F + o(l-e ) 


The main point to be made is that the assumption of independence is 
incorrect for the aircraft data since the sampling rate of 50 sanples/sec 
is 10 times as much as the Hyquist rate of the aircraft response. Hence 
the Kalman filter error variance will be larger than the nominal value p. 
The actual error variance 


E [(x - x ) 2 ] 

n n 


will depend upon the spectral density sliape of the noise process. V.’e shall 
calculate it for the sampled versions of Case 1 and Case 2. Thus 


satisfies: 


e = S A (1 - k ) e * (1-UrH (nA) + k n (n+lA) 
n+1 n r 
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Let 


R , = E C e 2 J. 
d n 


Case 1 
Let 

Y = B A 

Note that f > 1 if the sampling rate is s nailer than the Nyauist rate of 
2 

2B, and y < 1_ otherwise. In particular 

2 

Y = 0.05 

when the sampling rate is 10 times the Nyquist rate. We have 

R d = ( 1 ~ k) 2 T 2 k 2 o 2 Tan -1 ( 1 + T Tan Try ) 

1 - T 2 Try (1-T 2 ) 1 - T 


where 


T = ( 1 - k ) e 


-A 


To show the difference between the use of the continous-time model versus 
the discrete, we have plotted in Figure A-3, the ratio of to correspcn 

T 2 

error for the continous model for g" = o_ , as a function of y. Note 

2B 
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that as the sampling rate increases, the error on the discrete -model 
increases. The increase depends on the noise bandwidth B, increasing with 

B. 

Similar calculations carried out for Case 2 are plotted in Figure 

2 2 

A-4 where the continous-time model uses g = a . The curves are 

mfc 

similar. Note that again the substantial increase in error for the discrite- 
time model for y = .05, corresponding to the situation in our flight-test 


data. 
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